Data files from expt.(from Felix, 26/01/22): experimentally extracted DCS from one of our datasets, assuming a Up of 37.5eV and three different return energies of 1.2Up, 1.6Up, and 2.0Up, corresponding to scattering energies of 45eV, 60eV, and 75eV, respectively.
Data files for atomic scattering (from Felix, 26/01/22): ELESPA calculations for the C and H atoms and the three relevant energies. The columns are: scattering angle | real part of C scattering amplitude | imaginary part of C scattering amplitude | real part of H scattering amplitude | imaginary part of H scattering amplitude.
07/03/22 v1
Based on "EDCS $N_2$ Tests", 27/09/19, ePSproc_EDCS_N2_multiE_270919_dist.ipynb, see https://osf.io/nhdkf/
Basic read + plot for EDCS data using epsproc & Xarray.
(Note that only the dev brach of epsproc currently has the EDCS IO code required here.)
importsysimportosimportnumpyasnp# For module testing, include path to local module code here# modPath = r'/home/femtolab/github/ePSproc/'# sys.path.append(modPath)importepsprocasep
* sparse not found, sparse matrix forms not available.
* natsort not found, some sorting functions not available.
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available.
In [2]:
# Load data from modPath\data\filename# dataPath = r'/home/femtolab/python/remote/N2_tests/'# dataPath = '/home/paul/ePS/butadiene_2022/EDCSmultiE/orb15_BG/' # JakedataPath='/home/jovyan/projects/butadiene_2022/EDCSmultiE/orb15_BG/'# Docker hub version# Scan file(s), return list of Xarrays, one per file.dataSet=ep.readMatEle(fileBase=dataPath,recordType='EDCS')#, verbose = 3)
*** ePSproc readMatEle(): scanning files for EDCS segments.
*** Scanning dir
/home/jovyan/projects/butadiene_2022/EDCSmultiE/orb15_BG/
Found 1 .out file(s)
*** FileListSort
Prefix: /home/jovyan/projects/butadiene_2022/EDCSmultiE/orb15_BG/EDCSmultiE.orb15_BG_E10.0_10.0_150.0eV.inp.out
1 groups.
*** Reading ePS output file: /home/jovyan/projects/butadiene_2022/EDCSmultiE/orb15_BG/EDCSmultiE.orb15_BG_E10.0_10.0_150.0eV.inp.out
*** IO.fileParse() found 15 segments with
Start: EDCS - differential cross section program
End: ['+ Command', 'Time Now'].
Found 15 EDCS segments (sets of scattering results).
Processed 15 sets of EDCS file segments, (0 blank)
# Full map# ep.plot.hvPlotters.curvePlot(dataSet[0], kdims = 'Theta');dataSet[0].hvplot()#.hvPlotters.curvePlot(dataSet[0], kdims = 'Theta');
WARNING:param.RasterPlot02478: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
Out[4]:
In [5]:
# Line outs with widget dataSet[0].hvplot.line(width=1000)
Out[5]:
In [6]:
# Line outs with overlay# dataSet[0].pipe(np.log10).hvplot.line(width=1000, legend = 'left').overlay('E') # Legend postion not working?dataSet[0].hvplot.line(width=1000,xlim=[0,220]).overlay('E')# Add padding!
# Full map, log10 scaling# ep.plot.hvPlotters.curvePlot(dataSet[0], kdims = 'Theta');dataSet[0].pipe(np.log10).hvplot()
WARNING:param.RasterPlot07081: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
Out[7]:
In [8]:
# Line outs with widget dataSet[0].pipe(np.log10).hvplot.line(width=1000)
Out[8]:
In [9]:
# Line outs with overlay# ep.plot.hvPlotters.curvePlot(dataSet[0], kdims = 'Theta');# dataSet[0].pipe(np.log10).hvplot.line(width=1000, legend = 'left').overlay('E') # Legend postion not working?dataSet[0].pipe(np.log10).hvplot.line(width=1000,xlim=[0,220]).overlay('E')# Add padding!
Where $\sigma_{A}$ is the (incoherent) sum of the atomic DCSs.
Ito, Y. et al. (2017) ‘Extraction of geometrical structure of ethylene molecules by laser-induced electron diffraction combined with ab initio scattering calculations’, Physical Review A, 96(5), p. 053414. doi:10.1103/PhysRevA.96.053414.
Data files for atomic scattering (from Felix, 26/01/22): ELESPA calculations for the C and H atoms and the three relevant energies. The columns are: scattering angle | real part of C scattering amplitude | imaginary part of C scattering amplitude | real part of H scattering amplitude | imaginary part of H scattering amplitude.
# Set indexes# Grab Elist - should use parselinedigits regex here!Elist=[int(item.split('_')[-1].split('.')[0])foritemindataFileList[fileInds]]Elist# Set multindex for PD columnscols=pd.MultiIndex.from_product([Elist,['C','H']])cols
# Easiest to read with Numpy or Pandas, then push to Xr?importpandasaspd# This almost works, but seems to have issue with space vs. tab seperated files?# testList = [pd.read_csv(fIn, header=None, names=['Re','Im'], sep='\s+') for fIn in dataFileList[0:3]]testList=[pd.read_csv(fIn,header=None,index_col=[0],delim_whitespace=True)forfInindataFileList[fileInds]]# testXR = [pd.read_csv(fIn, header=None, names=['Re','Im'], sep='\s+').to_xarray() for fIn in dataFileList[fileInds]]# test# Sort-of works, but need to force/set dims first# E.g. cell 13 https://phockett.github.io/ePSdata/XeF2-preliminary/xe-xef2_comp_120221-dist.html# import xarray as xr# xrData = xr.concat(testXR, dim = 'E')# Try restacking from PDreformat=[pd.DataFrame([item[1]+item[2]*1j,item[3]+item[4]*1j]).TforitemintestList]# TODO tidy index, e.g. np.round(reformat[0].index,3)# [item.set_index(np.round(item.index,0), inplace=True) for item in reformat] # Tidy up index & round to nearest degree (to match ePS case, although interp might be better)# Set multindex# pd.MultiIndex.from_product(['C','H'],)pdData=pd.concat(reformat,axis=1)#, names = dataFileList[0:3])# pdData.columns = dataFileList[0:3] # Full paths at the moment - needs to tidy!# pdData.columns = [item.split('/')[-1].split('_')[-2] for item in dataFileList[fileInds]] # Quick and dirty enery settings# pdData.to_xarray() # OKpdData.index.rename('Theta',inplace=True)# pdData.rename(columns={'Eret1p2':45, 'Eret1p6':60, 'Eret2p0':75}, inplace=True) # Set energies from Up# pdData.rename(columns={0:'C',1:'H'}, inplace=True)pdData.columns=colspdData# Show table
Out[188]:
45
60
75
C
H
C
H
C
H
Theta
0.000000
1.255280e-08+1.383240e-08j
6.338500e-09+2.279280e-09j
1.366940e-08+1.305570e-08j
6.143260e-09+1.870990e-09j
1.442830e-08+1.254110e-08j
6.007380e-09+1.612820e-09j
0.000100
1.255280e-08+1.383240e-08j
6.338500e-09+2.279280e-09j
1.366940e-08+1.305570e-08j
6.143260e-09+1.870990e-09j
1.442830e-08+1.254110e-08j
6.007380e-09+1.612820e-09j
0.000125
1.255280e-08+1.383240e-08j
6.338500e-09+2.279280e-09j
1.366940e-08+1.305570e-08j
6.143260e-09+1.870990e-09j
1.442830e-08+1.254110e-08j
6.007380e-09+1.612820e-09j
0.000150
1.255280e-08+1.383240e-08j
6.338500e-09+2.279280e-09j
1.366940e-08+1.305570e-08j
6.143260e-09+1.870990e-09j
1.442830e-08+1.254110e-08j
6.007380e-09+1.612820e-09j
0.000175
1.255280e-08+1.383240e-08j
6.338500e-09+2.279280e-09j
1.366940e-08+1.305570e-08j
6.143260e-09+1.870990e-09j
1.442830e-08+1.254110e-08j
6.007380e-09+1.612820e-09j
...
...
...
...
...
...
...
178.000000
6.641090e-09-4.440460e-09j
1.113070e-10+9.907340e-10j
4.722510e-09-3.819200e-09j
1.398190e-10+7.038480e-10j
3.396300e-09-3.367850e-09j
1.451210e-10+5.384350e-10j
178.500000
6.643560e-09-4.440270e-09j
1.112030e-10+9.906440e-10j
4.724250e-09-3.818860e-09j
1.397450e-10+7.037790e-10j
3.397570e-09-3.367520e-09j
1.450650e-10+5.383810e-10j
179.000000
6.645320e-09-4.440140e-09j
1.111290e-10+9.905790e-10j
4.725480e-09-3.818630e-09j
1.396930e-10+7.037300e-10j
3.398470e-09-3.367290e-09j
1.450250e-10+5.383420e-10j
179.500000
6.646380e-09-4.440060e-09j
1.110850e-10+9.905400e-10j
4.726230e-09-3.818480e-09j
1.396620e-10+7.037010e-10j
3.399010e-09-3.367150e-09j
1.450010e-10+5.383190e-10j
180.000000
6.646730e-09-4.440040e-09j
1.110700e-10+9.905270e-10j
4.726480e-09-3.818440e-09j
1.396510e-10+7.036910e-10j
3.399190e-09-3.367100e-09j
1.449930e-10+5.383110e-10j
606 rows × 6 columns
In [199]:
# Need to concat with control here... this currently gives 6 das# atomicData = pdData.to_xarray() # Dataset# Use UNSTACK of course...atomicData=pdData.unstack().to_xarray().rename({'level_0':'E','level_1':'species'})atomicData
# atomicData.pipe(np.abs).hvplot.line(x='Theta').overlay('E').layout('species') # Abs (atomicData.real.hvplot.line(x='Theta').overlay('E').layout('species')+atomicData.imag.hvplot.line(x='Theta').overlay('E').layout('species')).cols(1)# Re + Imag
Out[221]:
TODO: these are amplitudes, so need to abs^2 and maybe renorm to match ePS DCS (angs^2).
In [224]:
# atomicData.pipe(np.abs).pipe(np.power,2).pipe(np.divide,(1E-10)**2).hvplot.line(x='Theta').overlay('E').layout('species') # Units x1000 compared to ePS?# atomicData.pipe(np.abs).pipe(np.power,2).pipe(np.divide,1E-16).hvplot.line(x='Theta').overlay('E').layout('species') # Units similar magnitudes to ePS?atomicData.pipe(np.power,2).pipe(np.divide,1E-16).pipe(np.abs).hvplot.line(x='Theta').overlay('E').layout('species')# Ordering might matter for coherent/incoherent results.# atomicData.pipe(np.abs).pipe(np.power,2).hvplot.line(x='Theta').overlay('E').layout('species') # Base units give ~1E-16, likely cm^2?# Q: scaling/normalisation here? Should expect similar magnitudes to ePS case?# If ELESPA in cm^2, get values 0-10 angs^2, maybe OK or a bit small?
# Do a quick and dirty interp for 5eV values - should be OK as not massive changes here, but ultimately should run higher-resolution ePS Eke scan.epsInterp=dataSet[0].sel(E=slice(35,85)).interp(E=np.arange(35,85,5),method='cubic')epsInterp.hvplot.line(x='Theta').overlay('E')+dataSet[0].sel(E=slice(35,85)).hvplot.line(x='Theta').overlay('E')
Out[136]:
In [298]:
# Run setPlotDefault for size etc.fromepsproc.plotimporthvPlottershvPlotters.setPlotDefaults(fSize=[1000,400])
In [299]:
# Calculate MCS # Square & rescale atomic data to angs^2# Sum over scatterers (incoherent)# Don't change Theta coords - Xarray should take care of that...? Should be OK as long as enough matching points (more angles in ELESPA outputs)sigmaAtomic=atomicData.pipe(np.power,2).pipe(np.divide,1E-16).pipe(np.abs)# Convert to XS in angs^2MCS=(epsInterp-sigmaAtomic.sum('species'))/sigmaAtomic.sum('species')MCS.hvplot.line(x='Theta').overlay('E')
Data files from expt.(from Felix, 26/01/22): experimentally extracted DCS from one of our datasets, assuming a Up of 37.5eV and three different return energies of 1.2Up, 1.6Up, and 2.0Up, corresponding to scattering energies of 45eV, 60eV, and 75eV, respectively.
# Easiest to read with Numpy or Pandas, then push to Xr?importpandasaspd# This almost works, but seems to have issue with space vs. tab seperated files?# testList = [pd.read_csv(fIn, header=None, names=['Re','Im'], sep='\s+') for fIn in dataFileList[0:3]]testList=[pd.read_csv(fIn,header=None,names=['Re','Im'],delim_whitespace=True)forfInindataFileList[fileInds]]testXR=[pd.read_csv(fIn,header=None,names=['Re','Im'],sep='\s+').to_xarray()forfInindataFileList[fileInds]]# test# Sort-of works, but need to force/set dims first# E.g. cell 13 https://phockett.github.io/ePSdata/XeF2-preliminary/xe-xef2_comp_120221-dist.html# import xarray as xr# xrData = xr.concat(testXR, dim = 'E')# Try restacking from PDreformat=[pd.DataFrame(item['Re']+item['Im']*1j)foritemintestList]# TODO tidy index, e.g. np.round(reformat[0].index,3)[item.set_index(np.round(item.index,0),inplace=True)foriteminreformat]# Tidy up index & round to nearest degree (to match ePS case, although interp might be better)pdData=pd.concat(reformat,axis=1)#, names = dataFileList[0:3])# pdData.columns = dataFileList[0:3] # Full paths at the moment - needs to tidy!pdData.columns=[item.split('/')[-1].split('_')[-2]foritemindataFileList[fileInds]]# Quick and dirty enery settings# pdData.to_xarray() # OKpdData.index.rename('Theta',inplace=True)pdData.rename(columns={'Eret1p2':45,'Eret1p6':60,'Eret2p0':75},inplace=True)# Set energies from UppdData# Show table
Out[302]:
45
60
75
Theta
1.0
0.000000+0.000000j
0.0000000+0.0000000j
0.000000+0.000000j
2.0
0.000000+0.000000j
0.0000000+0.0000000j
0.000000+0.000000j
3.0
0.000000+0.000000j
0.0000000+0.0000000j
0.000000+0.000000j
4.0
0.000000+0.000000j
0.0000000+0.0000000j
0.000000+0.000000j
5.0
0.000000+0.000000j
0.0000000+0.0000000j
0.000000+0.000000j
...
...
...
...
176.0
448.68140+55.12950j
168.412120+31.832830j
83.77806+20.52115j
177.0
452.79170+58.78483j
209.359720+40.423974j
53.27593+17.85868j
178.0
525.33100+73.77291j
226.295900+47.872022j
76.69515+23.49546j
179.0
639.52080+95.98866j
247.829700+56.554400j
201.85810+48.23733j
180.0
581.5425+115.8601j
236.015000+63.990280j
208.95260+62.50720j
179 rows × 3 columns
In [303]:
# pdData.to_xarray().real.hvplot() # Working for dataset for string col headers, but not ints?expData=pdData.to_xarray().to_array().rename({'variable':'E'})expData.pipe(np.abs).hvplot.line(x='Theta').overlay('E')# Abs
As above, but with scale factor for exp data. Should also be squared?
In [310]:
# Assume DCS - this looks interestingSF=1E-3# Basic order-of-magnitude scaling looks pretty good! # 4*7E-3 # SF should be E dependent? ALSO, most of the structure from the division, not subtraction.# SF = 1E-1 * 1/MCSexp.max('Theta') # E-dependent renormMCSexp=(SF*expData.pipe(np.abs)-sigmaAtomic.sum('species'))/sigmaAtomic.sum('species')# Squared? Less good?# SF = 6E-7# MCSexp = (SF*expData.pipe(np.power,2).pipe(np.abs) - sigmaAtomic.sum('species'))/sigmaAtomic.sum('species')# PlotMCSexp.hvplot.line(x='Theta').overlay('E')
